Bold Diagrammatic Monte Carlo Applied to Fermionized Frustrated Spins 
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We demonstrate, by considering the triangular lattice spin- 1/2 Heisenberg model, that Monte 
Carlo sampling of skeleton Feynman diagrams within the fermionization framework offers a universal 
first-principles tool for strongly correlated lattice quantum systems. We observe the fermionic sign 
blessing — cancellation of higher order diagrams leading to a finite convergence radius of the series. 
We calculate magnetic susceptibility of the triangular-lattice quantum antiferromagnet in the cor- 
related paramagnet regime and reveal surprising — numerically exact — microscopic correspondence 
with its classical counterpart at all accessible temperatures. Extrapolation of the observed relation 
to zero temperature suggests the absence of the magnetic order in the ground state. We critically 
examine implications of this unusual scenario. 
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The method of Bold Diagrammatic Monte Carlo 
(BDMC) [1] allows one to sample contributions from mil- 
lions of skeleton Feynman diagrams and extrapolate re- 
sults to the infinite diagram order, provided the series is 
convergent (or subject to re-summation beyond the con- 
vergence radius). Recent experimentally certified appli- 
cation of BDMC to unitary fermions down to the point 
of the superfluid transition [2] makes a strong case for 
BDMC as a generic method for dealing with correlated 
fermions described by Hamiltonians without small pa- 
rameters. One intriguing avenue to explore is to apply it 
to frustrated lattice spin systems, where, on one hand, 
standard Monte Carlo (MC) fails because of the sign 
problem [3], and, on the other hand, the system's Hamil- 
tonian can be always written in the fermionic represen- 
tation [4-6] which contains no large parameters — exactly 
what is needed for the anticipated convergence of BDMC 
with the diagram order. 

The BDMC approach is based on the sign blessing phe- 
nomenon, when, despite the factorial increase in the num- 
ber of diagrams with expansion order, the series features 
a finite convergence radius because of dramatic (sign al- 
ternation induced) compensation between the diagrams. 
With the finite convergence radius, the series can be 
summed either directly, or with re-summation techniques 
that can be potentially applied down to the critical tem- 
perature of the phase transition, if any. (At the criti- 
cal temperature thermodynamic functions become non- 
analytic, and the diagrammatic expansion involving ex- 
plicit symmetry breaking by the finite order parameter is 
necessary to treat the critical region and the phase with 
broken symmetry.) In the absence of the sign blessing, 
the re-summation protocols become questionable in view 
of the known mathematical theorems regarding asymp- 
totic series. At the moment, there is no theory allowing 
one to prove the existence of a finite convergence radius 
analytically. The absence of Dyson's collapse [7] in a 



given fermionic system is merely providing a hope that 
the corresponding diagrammatic series is not asymptotic 
and cannot be a priori taken as a sufficient condition for 
the sign blessing. Hence, the applicability of BDMC to 
a given system can be established only on the basis of a 
direct numerical evidence for series convergence and com- 
parison with either experiment or alternative controllable 
techniques, such as high-temperature series [8] and nu- 
merical linked cluster (NLC) expansions [9]. In this Let- 
ter, we report the first successful application of BDMC 
to fermionized quantum spin systems by simulating the 
canonical model of frustrated quantum magnetism — the 
triangular lattice antiferromagnetic spin- 1/2 Heisenberg 
model (TLHA). We demonstrate that BDMC for this 
frustrated magnet is indeed subject to the sign blessing 
phenomenon which allows us to obtain basic static and 
dynamic correlation functions with controllable (about 
one percent or better) accuracy. The agreement with 
extrapolated high-temperature expansions is excellent. 

In addition, we report a very surprising finding of ex- 
treme similarity between short-distance static spin corre- 
lations of the quantum and classical spin models, evalu- 
ated at different but uniquely related to each other tem- 
peratures. This numerically exact (within the error bars) 
quantum-to-classical correspondence holds at all accessi- 
ble to us temperatures T > 0.375 (here and below tem- 
perature is measures in the units of the exchange constant 
J). Specifically, the entire static correlation function of 
the quantum model at a given temperature T — having 
quite non-trivial pattern of sign- alternating spatial de- 
pendence and temperature evolution, thus forming a sys- 
tem's fingerprint — turns out to be equal, up to a global 
temperature dependent normalization factor, to its clas- 
sical counterpart at a certain temperature T c \ = T C \(T). 
Extrapolation of the obtained T C \(T) curve to T = limit 
results in a finite value of T c \(0) > 0, suggesting quantum- 
disordered ground state of the quantum model. 
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The Hamiltonian of the TLHA is given by 

H = JY J Si-Sj. (1) 

<i,3> 

Here Si is the spin- 1/2 operator on the z-th site of 
the triangular lattice and the sum is over the nearest 
neighbor pairs coupled by the positive exchange integral, 
J > 0. As found by Popov and Fedotov [4, 5], the grand 
canonical Gibbs distribution of the model (1) can be re- 
formulated identically in terms of purely fermionic oper- 
ators using 
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where is the second quantized operator annihilating 
a fermion with spin projection a,/3 = ±1 on site z, and 
a are the Pauli matrices. The representation (2) leads to 
a flat-band fermionic Hamiltonian, ifp, with two-body 
interactions and amenable to direct diagrammatic treat- 
ment. To eliminate statistical contributions from non- 
physical states having either zero or two fermions, Ref. 4 
introduced an imaginary chemical potential to H-p: 
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The added term commutes with the original Hamilto- 
nian and has no effect on properties of the physical sub- 
space {rii = 1} whatsoever. Moreover, the grand canon- 
ical partition functions and spin-spin correlation func- 
tions of the original spin model and its fermionic version 
are also identical because (i) physical and non-physical 
sites decouple in the trace and (ii) the trace over non- 
physical states yields identical zero on every site. As a 
result, one arrives at a rather standard Hamiltonian for 
fermions interacting through two-body terms. A complex 
value of the chemical potential, which can also be viewed 
as a peculiar shift of the fermonic Matsubara frequency 
uj n = 2?r(n + 1/2)T 2?r(n + 1/4)T, is a small price to 
pay for the luxury of having the diagrammatic technique. 

We perform BDMC simulations using the standard 
G 2 W- skeleton diagrammatic expansion of the fermionic 
model (l)-(3) in the real space-imaginary time repre- 
sentation [10], see also [11]. The first and most impor- 
tant question to answer is whether the sign blessing phe- 
nomenon indeed takes place. In Fig. 1 we show compari- 
son between the calculated answer for the static uniform 
magnetic susceptibility, Xu, and the NLC expansion re- 
sult [9] at T = 2. This temperature is low enough to 
ensure that we are in the regime of strong correlations 
because Xu is nearly a factor of two smaller than the 
free-spin answer Xu^ = 1/4T. On the other hand, this 
temperature is high enough to be sure that the high- 
temperature series can be described by Pade approxi- 
mants without significant systematic deviations from the 



exact answer [8, 9] (at slightly lower temperature the 
bare NLC series starts to diverge). We clearly see that 
the BDMC series converges to the correct result with ac- 
curacy of about three meaningful digits and there is no 
statistically significant change when more than a hun- 
dred thousand of 7-th order diagrams [12] are accounted 
for. The error bar for the 7-th order point is significantly 
increased due to factorial growth in computational com- 
plexity. The 4-th order result can be obtained after sev- 
eral hours of CPU time on a single processor. 

Interestingly enough, when temperature is lowered 
down to T = 1 which is significantly below the point 
where the bare NLC series start to diverge, see Fig. 2, 
the BDMC series continue to converge exponentially. 
This underlines the importance of performing simulations 
within the self-consistent skeleton formulation. 
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FIG. 1. (Color online) Uniform susceptibility calculated 
within the Gr 2 W-skeleton expansion as a function of the max- 
imum diagram order retained in the BDMC simulation (black 
dots) for T / J = 2. The result of the high-temperature expan- 
sion (with Pade approximant extrapolation) [8, 9] is shown by 
red square and horizontal line. 

In Fig. 2 we show results of the BDMC simulation 
performed at temperatures significantly below the mean- 
field transition temperature. We observe excellent agree- 
ment (within our error bars) with the Pade approximants 
used to extrapolate the high-temperature series data to 
lower temperature [8]. Within the current protocol of 
dealing with skeleton diagrams we were not able to go 
to lower temperature due to the development of near- 
singularity in the response function (and thus in the 
effective-interaction propagator) at the classical ordering 
wave- vector Q = (47r/3, 0), in units of inverse lattice con- 
stant. In future work we plan to apply pole-regularization 
schemes to overcome this technical problem. 

We now turn to the static susceptibility 
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FIG. 2. (Color online) Uniform susceptibility as a function 
of temperature (black dots) for the triangular Heisenberg an- 
tiferromagnet calculated within the BDMC approach. NLC 
expansion results based on triangles (labeled as 7T and 8T) 
and sites (labeled as 12S and 13S) [9] are shown along with 
two different Pade approximant extrapolations [8]. 



Here S z (r) is the Matsubara spin operator on the lat- 
tice site labeled by the integer index vector r. For 
simplicity of comparing susceptibility (4) with its clas- 
sical counterpart, we normalize it to unity at the ori- 
gin, x(r) — ^ x( r )/x(0)? doing the same with the classical 
Xci(r)- The latter is obtained by Metropolis simulation 
of the classical Heisenberg model (1) in which quantum 
spin operators are replaced with classical unit vectors n r . 
For every accessible temperature T we observe a perfect 
match (within the error bars) between quantum correla- 
tor x( r ) and its classical counterpart, for r ranging from 
1 to 5 (which includes 10 different sites), calculated at 
a certain temperature T C \(T). A typical example of the 
match is presented in Fig. 3. We note in passing that 
the equal-time correlation function, ( Sq(0) S*(0) ), while 
having qualitatively similar shape to that of (4), does 
not match the classical correlator Xci( r ) = ( n o n r)> es ~ 
pecially so for sites at which the sign of the correlation 
changes with temperature (such as sites 3 and 7 in Fig. 3 
for which the sign of x(r) changes from ferromagnetic at 
high T to antiferromagnetic one below T « 0.5). 

Mapping long-range correlations in quantum models 
onto the renormalized classical behavior is rather stan- 
dard on approach to the ordered state [13]. What we 
observe is fundamentally different: quantum-to-classical 
correspondence, or QCC, is valid in the intermediate 
temperature regime at all distances, including nearest- 
neighbor sites, and when the correlation length £ is still 
very short, of the order of the lattice spacing, £ ~ 1. It 
is worth noting that this short-distance correspondence 
is also very different from high-T quasi-classical wave 
regime of Ref. 13 which allows for the classical description 



at distances r ~ £ ^> 1. 

We find that QCC also takes place for the s = 1/2 
square lattice Heisenberg antiferromagnet where, thanks 
to the absence of sign problem for path-integral Monte 
Carlo, it can be verified with a better precision and over 
a much broader temperature range (down to the ground 
state in a finite-size system), 0.03125 <T/J< 1 and 1 < 
r < 5. These facts suggest that QCC is numerically exact 
and thus takes place at any temperature for both lattices, 
and, possibly, for other models of quantum magnetism. 

The quality of the matching procedure allows us to es- 
tablish the T C \(T) correspondence with accuracy of about 
one percent. In the right panel of Fig. 3 we plot the final 
result along with the asymptotic high-temperature rela- 
tion T c \ = (4/3) T reflecting the difference between the 
(S z ) 2 = 1/4 and ((n z ) 2 } = 1/3. An immediate conse- 
quence of observed QCC in Fig. 3 is that the entire q 
dependence of the static susceptibility x(q, uj n = 0) of 
the quantum model is given by the susceptibility of the 
classical model at temperature T c i(T), which is readily 
available from classical Monte Carlo simulations. 

Due to limited low-temperature range of the T C \(T) 
curve for the TLHA it is perhaps too early to make any 
definite conclusion regarding its extrapolation down to 
the T = limit. One possibility is that it smoothly ex- 
trapolates to a finite value T c i(0) = 0.28, implying that 
the ground state is some kind of a spin liquid. This possi- 
bility was discussed by Anderson [17] almost forty years 
ago but was subsequently rejected on the basis of nu- 
merous investigations which include exact diagonaliza- 
tion [18-20], Green's function MC [21], series expansion 
[22], density matrix renormalization group [23] stud- 
ies, as well as large-S (spin wave) [24-26], large-N [27] 
and functional renormalization group analysis [28]. Note, 
however, that the spin correlation length for the classi- 
cal model at T c \ « 0.28 is above 10 3 lattice periods [15] 
and thus simulations of small system sizes L ~ 10 would 
be severely affected by finite-size effects. The value of 
T c i(0) ~ 0.28 is surprisingly close (essentially within the 
error bars) to the temperature obtained by extrapolating 
transition temperatures for the q = 3 Potts transition 
in finite magnetic fields h to the h = limit [29, 30]. 
Large-scale MC simulations performed in zero magnetic 
field also identify T c \ = 0.285(5) as the critical point of 
the chiral transition [14-16]. However, the debate with 
regards to the existence of the chiral transition is not 
settled yet — an alternative scenario [31, 32] predicts a 
sharp crossover to a more standard non-linear sigma- 
model type behavior around T c \ = 0.28. 

The other possibility is that the QCC curve T c i(T) 
will cross over to the standard renormalized classical be- 
havior in the long wavelength limit and will arrive at 
T c \(0) = 0, implying the ordered quantum ground state. 
This is exactly what happens for the square lattice an- 
tiferomagnet, see inset in the right panel in Fig. 3. In 
fact, it is also known (and can be readily deduced from 
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FIG. 3. (Color online) Left and middle panels: Perfect match of the normalized quantum and classical responses at the 
corresponding temperatures, T and T C \(T). Points are ordered according to their distance from the origin, r, as is illustrated 
in the right top corner. The sign of the correlator is indicated explicitly next to the point. Right Panel: Mapping between the 
classical and quantum temperatures. The solid (continued as dashed) line is the fitting function, y — (4x 2 + Ax + B)/ (3x + C), 
with A — 0.462, B — 1.065, C — 3.825, which satisfies the asymptotic law y = 4x/3 expected in the high-temperature limit 
shown by the dashed-dotted line. Red arrow indicates the position of the chiral transition advocated in Refs. [14-16]. For 
comparison, the square lattice data (stars) are shown in the inset. 
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FIG. 4. (Color online) Blue (square) symbols: local static 
susceptibility x( r — 0,cj m = 0), multiplied by 4T, as a func- 
tion of T I J. Black circles show T-dependence of the (4 times) 
local spin correlation function x( r = 0,r = 13/2). 



the correspondence plot and Fig. 2 of Ref. 15) that in the 
TLHA the renormalized classical regime with large cor- 
relation length emerges only below temperature T « 0.25 
[33, 34], which is well below our lowest T = 0.375 data 
point. Clearly, more data at lower temperatures are re- 
quired in order to resolve this fascinating question. 

Normalization factor x(0) = x( r = 0>^m = 0) in Fig. 3 
is given by the local static susceptibility which of course 
is different for the classical and quantum system. For 
the classical Heisenberg model x(0) is simply 1/3, inde- 
pendent of temperature, while in the quantum system the 
local static susceptibility is T-dependent, as Fig. 4 shows. 



The same Figure also shows local spin correlation func- 
tion at r = 0/2, xir = 0,r = f3/2) = T^ m e^ m X (r = 
0,(jj m ). This too probes quantum fluctuations, i.e. con- 
tributions to the sum from terms with uj m ^ 0. As ex- 
pected, both curves deviate from unity with lowering of 
T, reflecting increasing role of quantum fluctuations. 

Perhaps the most striking feature of QCC is its pre- 
dictive power in the search for spin-liquid states. In- 
deed, if QCC is confirmed for a given model of quantum 
magnetism and the classical ground state is not ordered 
due to macroscopic degeneracy then the quantum ground 
state is not ordered as well, i.e. it is a spin liquid. (At 
this point we can only conjecture that this scenario takes 
place for the s = 1/2 Kagome antiferromagnet which re- 
cently has been shown to realize a long-sought spin 
liquid state [35].) Moreover, even if the classical ground 
state is ordered but the correspondence curve T C \(T) is 
such that T c i(0) ^ 0, the quantum ground state is still 
not ordered similarly to its finite-temperature classical 
counterpart. While the final outcome for the TLHA re- 
mains to be seen, our data convincingly show an unusual 
classical-to-quantum correspondence with regards to the 
static spin correlations. 
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